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To overcome the limitations of the traditional state-averaging approaches in excited state calcu- 
lations, where one solves for and represents all states between the ground state and excited state of 
interest, we have investigated a number of new excited state algorithms. Building on the work of 
van der Vorst and Sleijpen (SIAM J. Matrix Anal. Appl., 17, 401 (1996)), we have implemented 
Harmonic Davidson and State- Averaged Harmonic Davidson algorithms within the context of the 
Density Matrix Renormalization Group (DMRG). We have assessed their accuracy and stability of 
convergence in complete active space DMRG calculations on the low-lying excited states in the acenes 
ranging from naphthalene to pentacene. We find that both algorithms offer increased accuracy over 
the traditional State- Averaged Davidson approach, and in particular, the State- Averaged Harmonic 
Davidson algorithm offers an optimal combination of accuracy and stability in convergence. 



I. INTRODUCTION 

Many excited states possess complicated electronic 
structure which cannot be described by a single domi- 
nant electronic configuration. For such states, a reliable 
description requires a multireference quantum chemistry 
method. 

Recently, the Density Matrix Renormalization Group 
(DMRG) has emerged as a new tool for multireference 
quantum chemistry problems Q, 0, [], 4, 5, 6, 7]. When 
applied to bond-breaking, it achieves a balanced descrip- 
tion across potential energy curves due to its reference- 



conjugated polyenes and acenes [111. \12\. 

The DMRG ansatz can be written as a linear expansion 
in terms of many-body functions which are subsequently 
optimised with respect to internal non-linear degrees of 
freedom {R}, 



|*> = 5>ir|ir({R})> 



(1) 



Ir 



free nature [8l.l9l.ll0l|. Reduced-scaling DMRG algorithms 
have also been developed and applied to large multirefer- 
ence problems in quasi-one-dimensional systems such as 



Note that if we choose the expansion functions \lr) to be 
Slater determinants and the internal degrees of freedom 
{R} to be their constituent orbitals, the above ansatz de- 
scribes the Complete-Active-Space Self-Consistent-Field 

r~i 

(CASSCF) wavefunction [LJ. In the DMRG, the expan- 
sion functions are instead complicated many-body basis 
states and the non-linear degrees of freedom are renor- 
malisation matrices, which allows for a particularly com- 



pact and efficient expansion 



To obtain excited states in the DMRG we usually use 
the iterative Davidson algorithm to solve for eigenvec- 
tors = i>l r \lr) ranging from the ground-state to the 



excited state of interest 



15 



The non-linear parameters 



{R} for these states are subsequently optimised for a 
density matrix that is averaged over all the states 
State-averaging is necessary to improve the stability of 
the non-linear optimisation and to prevent root-flipping, 
which occurs when the approximate wavefunction leaves 
the convergence basin of the target excited state and en- 
ters that of a different excited state 
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The drawbacks of this conventional approach, which 
we shall refer to as the State-Averaged Davidson (SA-D) 
algorithm, become clear if one is interested in higher re- 
gions of the spectrum because it becomes infeasible, both 
in terms of computational cost and accuracy, to solve for 
and adequately represent all the lower-lying eigenvectors 
in the state-averaged DMRG basis. Consequently, it is 
desirable to explore alternative algorithms that directly 
yield individual or a few excited state wavefunctions at a 
time. Any such an algorithm should also retain the sta- 
bility of the SA-D algorithm during non-linear optimisa- 
tion, so as to be able to rapidly converge to the desired 
target excited state(s) without root-flipping. 

Iterative methods for linear algebra that work with 
shifted and inverted operators such as (ui — if) -1 have 
long been used in numerical analysis to obtain t he inte 
rior (i.e. excited state) eigenvalues of matrices 
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fication that used a shifted and inverted operator to di- 
rectly calculate harmonic Ritz approximations to excited 
eigenvalues and eigenvectors [23(. We shall refer to this 
variant as the Harmonic Davidson (HD) algorithm to 
distinguish it from the original algorithm above. Aside 
from a demonstration for the one-electron Kohn-Sham 



equation in Ref. 



24j , we are not aware of the application 



Sleijpen and van der Vorst proposed an efficient modi- 



of this technique elsewhere in quantum chemistry. 

The purpose of this work is to investigate the Har- 
monic Davidson algorithm as a means to directly tar- 
get individual excited states and regions of the spectrum 
within the DMRG. One area in which the current appli- 
cation to quantum chemistry differs from previous nu- 
merical applications is the presence of a subsequent non- 
linear optimisation step for the wavefunction. We inves- 
tigate how combining the Harmonic Davidson procedure 
with state-averaging over nearby states in the spectrum 
(State-Averaged Harmonic Davidson, or SA-HD) can be 
used to confer stability in this non-linear optimisation. 
While we have focused on the DMRG method here, our 
findings are relevant to excited state algorithms for other 
quantum chemistry methods whose ansatz contains both 
linear and non-linear parameters, such as in the CASSCF 
method. 

The structure of this paper is as follows. In Sec. UH we 
briefly review the DMRG method and the Davidson and 
Harmonic Davidson algorithms. In Sec. IIII1 we present 
DMRG calculations on the excited states of acenes from 
naphthalene to pentacene using both direct targeting 



with the Harmonic Davidson algorithm (in both state- 
averaged and non-state-averaged forms) as well as with 
the traditional (state-averaged) Davidson approach. We 
also compare our excited state spectrum with that ob- 
tained from Equation of Motion Coupled Cluster theory. 
We summarise our findings in Sec. IIVI 



II. THEORY 



A. DMRG 



The quantum chemistry DMRG algorithm used in this 
work has been described fully elsewhere 
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25J. As a 



detailed understanding is not necessary here, we shall 
restrict ourselves to only the essentials. As described 
above, the DMRG wavcfunction may be written in the 
form |T]). The DMRG sweep algorithm then provides 
an iterative method through which the many-body basis 
functions \l), \r) may be optimised with respect to a set of 
internal non-linear parameters R. For each orbital in the 
problem we can associate an R matrix, which describes a 
many-body renormalisation transformation involving the 
orbital (i.e. not simply an orbital rotation). In a sweep to 
optimize the |/) states (an analogous procedure holds for 
the |r) states), R matrices are determined from the M 
eigenvectors of the many-particle reduced density matrix 
with the largest eigenvalues. In the ground-state case, the 
density matrix that determines the \l) states is obtained 



by tracing out the \r) states from the wavefunction, viz 
Tw = ^ipiripi'r (2) 

r 

TwRi'm = li Rim, m=l,...,M (3) 

M is referred to as the size of the DMRG many-body 
basis, and as M increases, the DMRG wavefunction be- 
comes exact. For excited state calculations, it is usual 
to employ state-averaging to increase the stability of the 
non-linear optimisation. This consists of using an aver- 
aged reduced density matrix in eq. ([2]) 



(4) 



where typically we choose equal weights for all the states 
of interest. 

B. The Davidson Algorithm 

The Davidson algorithm provides an efficient iterative 
solver for the large number of linear coefficients in the 
expansion of the ground-state DMRG wavefunction ^ 



27] . |\&) is expressed in an auxiliary basis {rji} (gen- 
erated by the Davidson iterations) 

i 

h)=vl\lr) (6) 

The coefficients Cj are determined by left-projection with 

fail 



Y(r ]j \H-E\r h )c i = a 



(7) 



where E is the approximate expectation value 
(4>\H\tp)/ Each iteration of the Davidson al- 
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gorithm, generates a new basis function \rj) from the 
current trial solution \ip) via 



\r,) = {^g{H)-E)-\H-Em 



(8) 



which is then orthogonalised against and added to the 
subspace {r/i}. 

To obtain excited state eigenvectors, the simple gener- 
alization known as the block Davidson or Davidson-Liu 
algorithm 



28, 



29( is typically used. Here a residual vector 
is generated for each of the states from the ground-state 
up to the target excited state. Solution of the subspace 
eigenvalue equation ([7]) then yields successive approxi- 
mations to all eigenstates up to the excited state of in- 
terest. In the subsequent non-linear optimisation of the 
excited state in the DMRG algorithm, the eigenvectors 
obtained from the block Davidson algorithm (i.e. from 
the ground-state to the target eigenvector of interest) are 
all averaged together in the density matrix ([4]). We shall 
refer to this combined procedure as the State-Averaged 
Davidson, or SA-D algorithm. 

From the above, we see that the primary drawbacks of 
the traditional SA-D approach are (i) computational cost 

- we must solve for all the states between the ground-state 
and excited state of interest, and (ii) decreased accuracy 

- since a single set of non-linear parameters must now 
represent multiple states rather than a single state. 



C. The Harmonic Davidson algorithm 

To avoid the need to solve for the states below the 
excited state of interest as in the Davidson algorithm 
above, classic shift and invert methods map the target 
excited state of the Hamiltonian H onto the ground-state 
of a shifted and inverted operator f2 



n = h- 1 = (w - h)- 1 



(9) 



The Harmonic Davidson algorithm introduced by Slei- 
jpen and van der Vorst 



231 ] (see also Ref. [22| for a clear 
review) extends the Davidson algorithm to work with the 
operator £1 without the need to explicitly compute the 
operator inverse in eqn. (J9j) . Each iteration generates a 
basis {rji}, but now we expand the target excited state 
I*) in {H w r,i} 



(10) 



Left projection with {rjiH u \ yields a generalized eigen- 
value problem 

( Vj H u; \(H- 1 -E- 1 )\H a , m )c i =0 
52[{Vi\Hv,\rH)i ~ E-^njH^H^a = (11) 

i 

where E^ 1 is the current approximation to (uj — E)^ 1 . 
E u is known as a harmonic Ritz approximation to the 
corresponding eigenvalue of H^. From (jTTJ) , we see that 
solving the eigenvalue equation for H~ x in the subspace 
{H^Vi} i s equivalent to solving the eigenvalue equation 
for the non-inverted operator H u where the trial solu- 
tion is expanded in the basis and the coefficients 
are obtained by right projection using a different space 



{(jfo-Hu, |}. This suggests that subspace {77,} for eqn. (fTTj) 
can also be generated from the trial solution l^} through 
a Davidson-type iteration 

\ V ) = (diag(i^) - E' U )-\H U - ELM) (12) 

where here E' u refers to the expectation value 
(iplHuip) / (ip\ip) , which is distinct from E u appearing in 
eqn. (fTTJ) . 

While we could obtain the excited state eigenvalues 
and eigenvectors directly from the generalized eigenvalue 
problem pip , in practice it is numerically more stable 
to consider a slightly different form. By Schmidt orthog- 
onalization, we can construct an orthogonal decomposi- 
tion {fji\ of {H u T]i} such that (fjj H u \H u fji) = 8 0l . Re- 
expressing the eigenvalue problem in this basis gives 

J2((vAHu\Vi)-E- 1 S j i)c i = (13) 

i 

From eqn. (|13p we see that implementing the Har- 
monic Davidson algorithm requires only minor alter- 
ations to the traditional Davidson routine relating to the 
change in the subspace from {rji} to {fji}. In essence, 
there are only two additional steps: the subspace func- 
tions are first multiplied by H u , and second, they are 
Schmidt orthogonalized to yield 

In our later DMRG calculations, we will refer to the 
use of the above iterative procedure to solve for the lin- 
ear coefficients together with the non-linear optimisation 
of the many-body basis functions \r) without state- 
averaging, collectively, as the Harmonic Davidson algo- 
rithm (HD). 
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While the operator has the target excited state 
of interest as its ground-state eigenvector, stable conver- 
gence is not guaranteed in the non-linear optimisation. 
However, the formulation of the excited state problem 
as a ground-state minimization, albeit with a different 
operator Q, illustrates that root-flipping is really no dif- 
ferent from the poor convergence that may be found in 
difficult ground-state DMRG calculations. Consequently, 
the same procedures may be used to eliminate the con- 
vergence difficulty: either we can increase the size M of 
the DMRG basis or we can employ a state-average over 
the competing states. While we do not know a priori 
which states will cause convergence difficulties, it is rea- 
sonable to assume that they must lie energetically near 
our state of interest. We have thus implemented two 
types of State-Averaged Harmonic Davidson (SA-HD) 
algorithms. In the first (referred to as simply SA-HD) 
we average over the first n excited states of ft. These 
correspond to the n excited states that lie immediately 
above our target excited state in the spectrum of H. In 
the second, we average over the n states which lie clos- 
est (on either side) to the target excited state in the H 
spectrum. We refer to this variant algorithm as SA-HDa. 

The second variant (SA-HDa) is particularly suited to 
an alternative way of using the shift lu. Rather than 
choosing a shift to target a specific excited state, we can 
instead choose to find the excited states around a given 
shift. If stable convergence is not achieved, we simply 
then increase the number of states used in the SA-HDa 
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average until convergence is recovered. In this way, we 
can patch together the spectrum piece by piece by using 
successively higher shifts. 

III. APPLICATION TO ACENES 

We have investigated the low-lying states of the acene 
series ranging from naphthalene (2-acene) to pentacene 
(5-acene). In the following subsections, we describe 
the details of the computations (Sec. IIII Aj) . examine 
the excitation energies using the State- Averaged, Har- 
monic Davidson, and State-Averaged Harmonic David- 
son DMRG algorithms (Sec. IIII Bp . and finally use the 
(near-exact) DMRG results to assess the accuracy of the 
excitation spectrum obtained from Equation-of-Motion 
Coupled Cluster theory (EOM-CC) (Sec. HITC|) . 

A. Computational Details 



We used a model geometry for the acenes with C% n 
symmetry. The C-H bond lengths were 1.090 A. Along 
the legs of the acene ladder, the alternate C-C bond 
lengths were 1.410 A and 1.405 A, respectively. Along 
the rungs of the acene ladder, the C-C bond length was 
1.465 A. An example geometry for naphthalene is shown 
in Fig. [T] 

All calculations used the Slater-Type-Orbitals fitted 
to 3 Gaussians minimal basis set (STO-3G), consist 
ing of 2slp functions on C and Is functions on H 




We obtained the atomic orbital integrals and Restricted 



FIG. 1: Naphthalene model geometry. 

Hartree-Fock (RHF) orbitals from the PSI3.2 package 
3 11 ] . The RHF energies are given in Table U For the 
excited state calculations, we used a 7r-active space con- 
sisting of one p z orbital per carbon i.e. n-acene would 
have a (4n + 2, An + 2) active space. In the DMRG cal- 
culations, we further symmetrically orthonormalized the 
p z orbitals with respect to the overlap S. This gave a 
local orthonormal basis which yields faster convergence 
in the DMRG calculations. The remaining non-active or- 
bitals from the RHF calculations were kept frozen in all 
calculations. 

We calculated excitation energies with the State- 
Averaged Davidson (SA-D), Harmonic Davidson (HD), 
and State- Averaged Harmonic Davidson (SA-HD) algo- 
rithms described in Sec. |TTJ Our calculations used the lo- 
cal quadratic-scaling DMRG algorithm described in Ref. 

We employed a screening threshold of 10 -8 Hartrees 
{Eh) with no spatial symmetry. The ordering of the or- 
bitals for anthracene is shown in Fig. [5] and the other 
acenes were ordered similarly. In all of our sweeps, we 
added a small amount of random noise (10~ 6 — 10~ 8 ) to 
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TABLE I: RHF, CCSD, and DMRG(500) total energies of the acenes. All energies are in hartrees. 



Molecule 


HiRHF UUoJJ DMnAj^oUUJ 




-378.66597 -378.85130 -378.85360 


C14H10 


-529.44420 -529.70634 -529.71032 


C18H12 


-680.21823 -680.56059 -680.56538 


C22H14 


-830.99045 -831.41614 -831.42016 



the density matrix so that we would not lose important 
quantum numbers 



25 



32 1 . In the current algorithm it 
is difficult to converge DMRG energies beyond the in- 
trinsic accuracy associated with the finite number M of 
DMRG basis states. Thus DMRG energies were con- 
verged to within 1 milliHartree (mEh) (M = 50), 0.5 
mE h (M = 100), 0.5 mE h (M = 250), or 0.1 mE h 
(M = 500), respectively. We note that our largest M 
DMRG excitation energies are essentially exact (within 
the one-particle basis) to all reported digits. This is pos- 
sible for the large active spaces used here because of the 
compact parametrisation afforded by the DMRG wave- 
function. 





FIG. 2: The orbital ordering used for anthracene. 

In the HD and SA-HD calculations, the shift u> for 
a specific root was obtained as follows. To begin, we 
guessed an initial shift (typically based on our previous 
SA calculations) . In the case where the shift was too low 



or too high, the next guess for ui was obtained from the 
DMRG (block) iteration, where an undesired state first 
appeared as the ground state of the Harmonic Davidson 
procedure. The shift lo was then taken to lie on the cor- 
rect side of the desired state in this iteration. In this 
simple manner, we found that we could obtain a suitable 
shift for a given root with at most two to three guesses. 

To determine the symmetries of the excitations in 
the DMRG calculations we used the following method. 
Firstly, spin symmetries were obtained from the expec- 
tation value of S 2 . To obtain the spatial symmetries, we 
first assumed that the ground-state 'J'o was of A\ symme- 
try (as in experiment). For the excited states, we exam- 
ined "dipole" type matrix elements (^ol^o + n o ~ n i ~ 
niY&i) (essentially a dipole transition element along the 
short-axis of the acene; and 1 refer to atom labels in 
Fig. [21) For singlet excited states a non- vanishing dipole 
then implied B2 symmetry, while a vanishing dipole im- 
plied Ai symmetry. For the triplet excited states, all 
such matrix elements vanish. However, we could still de- 
termine the spatial symmetry through the expectation 
value (^ol^-o — n i smce n o ~ n i does not preserve 
spin symmetry and creates a residual expectation value 



s 



from which one can determine the spatial symmetry of 
the excited state. 

To obtain the orbital character of the excitations, 
we calculated transition onc-particlc density matrices 
{^O where denotes the ith excited state and 
identified the largest matrix elements. 

We further calculated the excitation spectrum (in the 
same 7r-active space as the DMRG calculations) with 
the Equation-Of-Motion Coupled Cluster Singles and 

3J 1 usir 



Doubles method (EOM-CCSD) 



package 
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using the dalton 



B. Comparison of Excited-state algorithms for 
DMRG by SA, HD, and SA-HD 

The ground state DMRG energies for the acenes are 
given in Table [J Tables Hfl UHl HY) and E] contain the 
first seven 7r — 7r* excitation energies for each acene, 
while Fig. [3] displays them in graphical form. Under 
C2 V symmetry, the only two possible representations of 
the 7r — 7r* excited states are 1,3 A\ and 1,3 Ba. Experi- 
mentally, there are three well-documented singlet bands 
that appear in the visible spectrum 
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The a- 



band and /3-band correspond to a polarization along the 
long axis and the p-band corresponds to a transverse 
polarization. We observed the a-transition as the low- 
est singlet excitation in each acene. Neither the p-band 
nor the /3-band appeared within the first eight states of 
each acene. Instead, for the case of naphthalene, the 
p-band emerged at 8.42 eV (state 19). The p-band nor- 



mally appears lower in the spectrum, but the absence 
of dynamic a — t: correlations is responsible for its artifi- 
cially high excitation energy here. This is consistent with 
previous studies of acenes using Complete- Active-Space 
Self-Consistent-Field (CASSCF) and Complete-Active- 
Space Moller-Plesset second order perturbation theory 
(CASMP2) theory 
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38 



39j. Triplet excitations are 



somewhat harder to measure experimentally. We observe 
that the triplet excitation energies decrease in energy 
more rapidly with system size than the singlet excita- 
tions. Thus while in naphthalene and anthracene there 
is one triplet level between the first two singlet excita- 
tions, in naphthacene and pentacene there are two. 

Comparing the accuracies of the SA-D, HD, and SA- 
HD calculations we observe that as expected, (other than 
by the size of the DMRG basis M), the accuracy in the 
excitation energies is determined primarily by the num- 
ber of eigenvectors in the state-average. Consequently 
the traditional SA-D algorithm yielded the lowest accu- 
racy (as it averages over all states between the ground 
state and excited state of interest) while the HD calcula- 
tions were correspondingly the most accurate since they 
targeted a single state at a time. The accuracy of the 
SA-HD calculations lay somewhere in between depend- 
ing on the number of states used in the average. In all 
cases, the differences between the various algorithms was 
most marked for the smaller sizes M of the DMRG basis, 
as for larger M all the wavefunctions become essentially 
exact. We would expect the differences to become more 
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TABLE II: DMRG excitation energies for naphthalene (CioHg) obtained with the SA-D, HD and SA-HD algorithms. All 
energies are in eV. State refers to the ground-state, and SA[m-n] refers to a state-average over all states from the mth to 
nth excited state. Numbers in parentheses give the number of DMRG states M. The "Exact (HD(500))" numbers are the 
(near-exact) excitation energies, while other entries give the errors from this result. The "Excitation" row gives the character 
of the excitation where 1 denotes HOMO, 2 denotes HOMO-1, 1' denotes LUMO, 2' denotes LUMO+1 and so on. The last 
column gives the mean improvement in the excitation energy over the SA [0-7] D result with the same M. n.c. denotes no 
convergence. 
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sa fo-7i n ci nn"i 

uir 1 U ( J i-J 1 -L UKJ J 
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0.02 
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0.02 




SA [0-3] D (250) 


0.00 


0.00 


0.00 


0.00 










0.00 


SA [3-7] HD (250) 
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0.00 


0.00 



pronounced in larger systems, where we are unable to use 
a sufficiently large M to reach exactness. 

Regarding the stabilities of the various algorithms, we 



found that there were no difficulties in converging the 
DMRG sweeps to the correct states with the SA-D al- 
gorithm. The HD algorithm on the other hand exhib- 
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TABLE III: DMRG excitation energies for anthracene (C14H10). Refer to table ITT1 for details. 
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bA [U-7J JJ (1U0) 


0.12 


0.12 


0.15 


0.15 


0.25 


0.23 


0.18 


0.19 




oA [U-oJ JJ [J-UUJ 


n n7 
U.U ( 


a n7 
U.U ( 


a 1 a 
U.1U 


0.09 










A AK 
U.UO 


SA [3-7] HD (100) 








0.15 


0.20 


0.20 


0.16 


0.28 


0.00 


tj"pi /inn'i 
xlJJ (^J-UUJ 


U.Ul 


a aq 
U.Uo 


U.Uo 


0.04 


n.c. 


n.c. 


n.c. 


n.c. 


A 1 A 
U.1U 


oi\ [D-OJ HU ^lUUJ 












n 10 
U. lo 


U.ZZ 




n hq 
u.uo 


oA [D- / J ilU ^lUUJ 














U. 14 


U. 1Z 


u.uy 


SA [0-7] D (250) 


0.01 


0.01 


0.01 


0.02 


0.03 


0.02 


0.02 


0.02 




oA [U-o| JJ l^zoUJ 


n nn 
U.UU 


n nn 
U.UU 


U.U1 


0.01 










A Al 

U.Ul 


bA [o-/J 11JJ [ZOU) 








0.02 


0.03 


a no 
U.Uz 


A AO 

U.Uz 


A AO 

u.uy 


A Al 

— U.Ul 




0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.02 


SA [0-7] D (500) 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 




SA [0-3] D (500) 


0.00 


0.00 


0.00 


0.00 










0.00 


SA [3-7] HD (500) 








0.00 


0.00 


0.00 


0.00 


0.00 


0.00 



ited the expected convergence difficulties characteristic of 
root-flipping for certain higher excited states. As previ- 
ously discussed, the stability of the HD algorithm would 
increase with the size of the DMRG many-body basis 
Al. In naphthacene, we required M > 250 to converge 



states 5-7 with the HD algorithm, while in pentacene, 
we required M — 500 to converge states 4-7. While the 
HD algorithm exhibited root-flipping, it was ameliorated 
with respect to simple eigenvector following (defined as 
following the n th eigenvector in the block Davidson al- 
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TABLE IV: DMRG excitation energies for naphthacene (C18H12). Refer to table|II]for details. 



Method 








State 








Mean 




1 Ai 


1 -D2 




2Mi 


l 3 Ai 


ol A 

3 Ai 


O -D2 


o3 A 


Improvement 






1 _* 1' 

1 — c 1 


O — > 1 


2 -> 1' 


2-> 1' 


. o 7 

— > z 


P. _v 1 ' 

— > 1 


A v. 1 1 

4 — > 1 




Excitation 




O — > O 


1 V 

1 — » O 


1 -> 2' 


1 -> 2' 


z — > 


1 — > 


1 _^ A' 




li<xact^tiJJ^OUU ) ) 


n nn 
U.UU 


1 PCO 
l.OZ 


z.yo 


3.27 


3.50 


Q HQ 

o.yo 


A AO 

4.Uz 


A OQ 
4.ZO 




o A rn 71 T\ ^p;n\ 

oA [U- r j jj (oUj 


n 71 
U. ( 1 


U.ol 


l.U / 


1.33 


1.75 


1.40 


1 p; i 

l.ol 


i 7n 




OA rn Ol T\ /rrn\ 

bA [U-oJ JJ (oUj 


n Kn 
U.oU 


r\ ao 
U.4o 


U.OO 


0.92 










U.34 


oA [ZWJ xlJJ (^OUJ 






a on 
U.zU 


1.04 


1.25 


1 on 


1 Q /I 

l.o4 


1 OA 


n qq 
U.00 


SA [3-7] HD (50) 








n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


0.00 


rlJJ (OUJ 


n on 
U.zU 


n 07 
U.Z ( 


n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


A CO 
U.OZ 


oA [l-^j -H-JJ 




1.4o 


U.yo 












n or 
— U.zo 


OA To Ql U"Ti fp;n\ 
bA [z-oj xlJJ [oUJ 






n ^7 
U.O / 


0.84 










n a k 
U.4o 


OA rn 71 "Pi / 1 nr^ 
bA [U- rj JJ [tUUJ 


n 07 
U.z ( 


n oq 
U.zo 


U.oz 


0.40 


0.58 


A A 7 
U.4/ 


A A 1 

U.41 


n rn 
U.OU 




SA [0-3] D (100) 


0.12 


0.14 


0.16 


0.20 










0.16 


bA rlJJ [luUJ 






n qq 
U.00 


0.40 


0.57 


A /I K 
U.40 


A on 
U.oy 


U.04 


n aa 
U.UU 


bA [0- / J ilU ^1UUJ 








n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


n nn 

U.UU 


xliJ ^lUUJ 


u.uo 


u.uo 


n ns 

U.Uo 


0.11 


n.c. 


n.c. 


n.c. 


n.c. 


u.zo 


SA [0-7] D (250) 


0.03 


0.04 


0.04 


0.06 


0.08 


0.08 


0.06 


0.07 




OA TO Ql Tt /'OP.n^ 
oA [U-o| JJ [ZoU) 


U.U1 


n m 
U.U1 


n no 
U.Uz 


0.03 










n aq 
U.Uo 


OA TQ 7l UTi fOKf\\ 

bA [o-/J rlJJ [zOUJ 








0.05 


0.07 


n fi7 
U.U < 


n nfi 
U.UD 


n nn 

u.uy 


n aa 
U.UU 


HD OW) 

±±±y y^iOKJ j 


0.00 


0.00 


0.00 


n.c. 


n.c. 


0.02 


0.01 


0.01 


0.05 


SA [0-7] D (500) 


0.00 


0.00 


0.00 


0.01 


0.01 


0.01 


0.01 


0.01 




SA [0-3] D (500) 


0.00 


0.00 


0.00 


0.00 










0.00 


SA [3-7] HD (500) 








0.01 


0.01 


0.01 


0.01 


0.01 


0.00 



gorithm in successive DMRG iterations) because of the 
use of the shift uo. For example, with M — 100, the 
third excited state of naphthalene could not be converged 
with simple eigenvector following, but could be converged 
without difficulty using the HD algorithm. 



Including a sufficient number of states in the SA-HD 
algorithm restored the stability of the convergence. Cer- 
tain "competing" states were particularly important for 
the state average, especially for smaller M. For all the 
acenes, the second and third excited states were examples 
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TABLE V: DMRG excitation energies for pentacene (C22H14). Refer to table|IT]for details. 



Method 








State 








Mean 




1 Ai 


1 3 B 2 


2 3 B 2 


2Mi 


l 3 Ai 




3 3 B 2 


2 3 Ai 


Improvement 






1 — > 1 


2^1 


3 -> 1' 


3 -> 1' 


3^2' 


4 — » 1 


1 rf 
1^5 




Excitation 




z — > z 


1 — > z 


1 -> 3' 


1 -> 3' 


z — > 


1 A' 

1 — > 4 

. o' 

Z — r Z 


— r 1 




TT-vojrf fVfT\f KnriW 
ii/xact^J4JJ^0UU ) ) 


n nn 
U.UU 


1.10 


QO 

z.oy 


3.10 


3.15 


0.0U 


Q /I Q 
o.4o 


O.OO 




OA rn 71 "P\ / c n "\ 

bA [U-rJ JJ (oU) 


1 in 


1.00 


1 7H 


1.86 


2.31 


on 
z.29 


Z.ZD 


z.oo 




OA rn ol T~\ / trn\ 

bA [U-oJ JJ (5(J) 


0.72 


0.72 


0.98 


1.24 










0.66 


SA [2-7] HD (50) 






1.23 


1.62 


1.88 


1.75 


1.97 


1.95 


0.41 


oA [0- f J xiJJ [OU ) 








n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


n nn 
U.UU 


un /en'* 
J4JJ l^oUJ 


U.zy 


U.4U 


11. c. 


n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


U.yo 


bA [z-o| xiJJ [OUJ 






n 

U.Oo 


1.21 










U.00 


OA rn 71 "Pi / 1 nn"\ 
bA [U- / J JJ 


U.44 


n /id 
U.4o 


n co 
U.OZ 


0.70 


0.87 


n qo 
U.oU 


n co 
U.oz 


n 07 
U.o ( 




SA [0-3] D (100) 


0.31 


0.33 


0.34 


0.38 










0.20 


OA TO 71 "PTT\ /"inn 1 * 

bA riJJ [1UUJ 






U.4/ 


0.56 


0.75 


U.O / 


n 

U.OO 


n en 
U.Oy 


n 1 c 

U.10 


bA [0- ( J ilU ^1UUJ 








n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


n nn 

U.UU 


xliJ ^lUUJ 


n ri/i 

U.U4 


n nn 




n.c. 


n.c. 


n.c. 


n.c. 


n.c. 


n qq 
u.oy 


SA [0-7] D (250) 


0.06 


0.08 


0.10 


0.12 


0.15 


0.16 


0.12 


0.15 




OA rn Ql "n foKr\\ 
bA [U-o| JJ l^zoUJ 


n no 
U.U2 


n no 
U.Uz 


U.Uo 


0.04 










n nfi 

U.UO 


bA [oWJ riJJ [zOUJ 








0.09 


0.10 


n 1 
U.lz 


n 1 n 
U.1U 


n 1 
U.lz 


n no 
U.Uo 


HD OW) 

±±±y y^iOKJ j 


0.00 


0.01 


0.01 


0.02 


n.c. 








0.08 


SA [0-7] D (500) 


0.00 


0.01 


0.01 


0.01 


0.02 


0.02 


0.02 


0.02 




SA [0-3] D (500) 


0.00 


0.00 


0.00 


0.00 










0.01 


SA [3-7] HD (500) 








0.01 


0.01 


0.02 


0.02 


0.02 


0.00 



of such states. Thus while the state averages SA[2-3] HD 
and SA[2-7] HD converged without difficulty, calculations 
using SA[3-7] HD did not, at least for smaller M . 

As mentioned previously, rather than choosing a shift 
to target specific excited states, we could take the differ- 



ent approach of trying to find the excited states around 
the frequency of a given shift lo. In this way, we could 
piece together a complete spectrum by performing, say, 
SA-HD or SA-HDa calculations with successively higher 
shifts. To demonstrate this, we computed the excita- 
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TABLE VI: DMRG excitation energies for the higher excited states of naphthalene (CioHg). Refer to table [H] for details. 



Method 






State 






Mean 




2 3 Ai 


3Mi 


4 s B 2 


3 3 Ai 


1 1 B 2 




Improvement 




A 1 1 

4^1 


3^2' 


4^2' 


2^3' 


1^3 


4^1 




Excitation 


1 — ► 4 


2 — > 3' 


2^4' 


3^2' 


3 — » 1 

4 — > L 


1 — » 4 




Exact(HD(500)J 


5.51 


5.87 


6.28 


6.48 


6.84 


6.84 




of\- [U-11J 1 ) 


U.Zy 


U.Zl 


0.19 


0.45 


n qo 


U.4o 




jA [0-1 1J \\\ ) < I 


n oq 


n on 
u.zu 


0.17 


0.46 


u.oo 


n a q 


n nn 

U.UU 


xlJJ {0\J) 


n.c 


n oq 
U.U8 


n.c 


n.c 


11. c 


n.c 


n i q 
U.lo 


SA [0-11] D (100) 


0.03 


0.03 


0.03 


0.03 


0.05 


0.06 




oA [o-llj HJJa l^lUUj 


U.Uo 


n no 
U.Uz 


0.03 


0.03 


n n^ 


n n^ 
U.Uo 


n nn 
U.UU 


rliJ (^iUUJ 


U.U1 


n ni 
U.U1 


0.01 


0.01 


n no 
U.U2 


n.c 


n no 


bA [iU-11] xiJJa (^iUUJ 










n n<i 
U.U4 


n n^i 
U.U4 


n no 
U.Uz 


SA [0-11] D (250) 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 




oA [0-11J HUa yZo\J) 


U.UU 


n no 
U.UU 


0.00 


0.00 


n nn 
U.UU 


n nn 
U.UU 


n nn 
U.UU 


J. ±J-/ 1 ZjOKJ j 


0.00 


0.00 


0.00 


0.00 


0.00 




0.00 


SA [10-11] HDa (250) 










0.00 


0.00 


0.00 


SA [0-11] D (500) 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 




SA [6-11] HDa (500) 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 



tion energies for states 6-11 for naphthalene using the 
SA-HDa algorithm with a shift chosen slightly above the 
state 7 excitation energy as estimated from the previous 
SA-HD [4-7] calculation. These are shown in table IVll 



C. Comparison of DMRG and EOM-CC excitation 
energies in the acenes 

The ground state EOM-CCSD energies for the acenes 
are summarized in Table |TJ We used our near-exact 
DMRG(500) excitation energies to examine the accuracy 
of the EOM-CC method in acenes. The EOM-CCSD and 
the DMRG symmetries and excitation energies are shown 
in Fig. OH For the larger acenes, the EOM-CCSD excited 



states are in a qualitatively different order as compared 
to DMRG. Similarly, EOM-CCSD erroneously predicts a 
very small singlet-triplet gap for the longer acenes. This 
points to the necessity of including relatively high order 
correlation effects to accurately describe excitations in 
the acenes. 

IV. CONCLUSIONS 

To overcome the computational and accuracy limita- 
tions of the traditional State-Averaged Davidson Algo- 
rithm, which requires both solving for and representing 
all states between the ground state and excited state of 
interest, we have investigated a number of new excited 
state algorithms within the context of the Density Ma- 
trix Renormalization Group (DMRG). In the Harmonic 
Davidson (HD) algorithm, using a shifted and inverted 
operator enabled us to directly solve for the excited state 
of interest. In the State- Averaged Harmonic Davidson 
(SA-HD) algorithm, we combined the HD method with 
an average over nearby excited states, to confer greater 
stability and overcome problems of root-flipping in the 
non-linear optimisation of the wavefunction. 

To assess the accuracy, stability, and computational 
cost of these new methods we calculated the low-lying 
excited states in the acenes ranging from naphthalene to 
pentacene. We found that as expected, in addition to 
the size of the DMRG basis M used, the accuracy was 
primarily determined by the number of states used in 
the state average. Thus the State-Averaged Davidson 
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approach gave the least accuracy, the Harmonic David- 
son algorithm, the highest, and the State- Averaged Har- 
monic Davidson lay in between depending on how many 
nearby states were included. The State-Averaged Har- 
monic Davidson algorithm converged smoothly without 
root-flipping so long as nearby "competing" states were 
included in the average. 

We also argued that through the shift ui in the Har- 
monic Davidson algorithms we could piece together a 
complete excitation spectrum by targeting different re- 
gions with successively higher shifts. This we demon- 
strated by calculating some higher lying excited states in 
naphthalene. 

Within the basis used, our DMRG excitation ener- 
gies are near-exact and we have used them to assess 
the accuracy of the EOM-CCSD method in the acenes. 
We found that the EOM-CCSD excitation spectrum was 
qualitatively different from that of the DMRG for the 
larger acenes, which demonstrates the necessity of in- 
cluding higher-order correlations to properly describe the 
electronic spectrum of conjugated quasi-one-dimensional 
molecules. 

Finally, we observe that the Harmonic Davidson algo- 
rithms studied here are quite general methods and are not 
limited to the Density Matrix Renormalisation Group. 
Thus they may be useful also to target excited states in 
other multi-reference theories, such as Complete Active 
Space Self-Consistent-Field theory. 
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FIG. 3: Comparison of DMRG and EOM-CCSD excitation 
energies for acenes. All energies in eV. 



